#include <fstream>
#include <iostream>
#include <iomanip>
#include "TFile.h"
#include "TTree.h"
#include "TLine.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TPaveText.h"
#include "TNtuple.h"
#include "TGraph.h"
#include "TString.h"
#include "TArc.h"
#include "TEllipse.h"
#include "TClonesArray.h"
#include "../Objects/ObjectClass.hh"

using namespace std;

void gamma2(){

  int ret = gSystem->Load( "../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }

  TFile *fread=new TFile("/scratch/shinichi/OutputFiles/gamma_5000.root");
  if(!fread){ cout<<"No Such File"<<endl; return;}

  TTree *tree=(TTree*)fread->Get("tree");
  
  TClonesArray *JPMTa=new TClonesArray("PMT");
  TClonesArray *MPMT1a=new TClonesArray("PMT");
  TClonesArray *MPMT2a=new TClonesArray("PMT");
  TClonesArray *MPMT3a=new TClonesArray("PMT");
  TClonesArray *MPMT4a=new TClonesArray("PMT");
  TClonesArray *Particlea=new TClonesArray("Particle");
  TClonesArray *Neutrona=new TClonesArray("Particle");
  TClonesArray *Argona=new TClonesArray("Particle");
  TClonesArray *Gammaa=new TClonesArray("Particle");
  int nPhoton;
  int nElectron;

  tree->GetBranch("JPMTb")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT1b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT2b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT3b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT4b")->SetAutoDelete(kFALSE);
  tree->GetBranch("Particleb")->SetAutoDelete(kFALSE);
  tree->SetBranchAddress("JPMTb",&JPMTa);
  tree->SetBranchAddress("MPMT1b",&MPMT1a);
  tree->SetBranchAddress("MPMT2b",&MPMT2a);
  tree->SetBranchAddress("MPMT3b",&MPMT3a);
  tree->SetBranchAddress("MPMT4b",&MPMT4a);
  tree->SetBranchAddress("Gammab",&Gammaa);
  tree->SetBranchAddress("nPhoton",&nPhoton);
  tree->SetBranchAddress("nElectron",&nElectron);

  const int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;

  TFile *fwrite=new TFile("out_5000.root","recreate");

  TNtuple *nt[nev];
  TH1F *hNPhoton=new TH1F("hPhoton","hPhoton",1000,0,100000);
  TH1F *hNElectron=new TH1F("hElectron","hElectron",500,0,50000);

  int cFid=0;
  int cFidInteraction=0;

  for(int iev=0;iev<nev;iev++){
    JPMTa->Clear();
    MPMT1a->Clear();
    MPMT2a->Clear();
    MPMT3a->Clear();
    MPMT4a->Clear();
    Gammaa->Clear();

    nt[iev]=new TNtuple(Form("nt_%d",iev),Form("nt_%d",iev),"itrack:istep:isFid:ProcessBit:xpre:ypre:zpre:xx:yy:zz:edep:ekin");

    tree->GetEntry(iev);

    bool IsFid = false;

    //----- Gamma
    int igamma = 0;
    int trackid=-1;
    int istep = 0;
    cFidInteraction = 0;
    for(int i=0;i<Gammaa->GetEntriesFast();i++){
      Particle *gamma = (Particle*)Gammaa->At(i);
      TString processname = gamma->processname; 
      TString volumename = gamma->volumename;
      int trackID = gamma->TrackID;
      float xpre = gamma->xpre;
      float ypre = gamma->ypre;
      float zpre = gamma->zpre;
      float xpos = gamma->xpos;
      float ypos = gamma->ypos;
      float zpos = gamma->zpos;
      float edep = gamma->edep;
      float ekin = gamma->ekin;

      if(trackid!=trackID){
	trackid = trackID;
	istep = 0;
	igamma++;
      }

      if( processname!="Transportation"&&volumename=="pFid"){ 
	IsFid=true;
	cFidInteraction++;
      }

      int ProcessBit = 0;
      int isFid   = 0;
      if( processname=="Transportation" ) ProcessBit = 1;     
      if( processname=="compt" ) ProcessBit = 2;     
      if( processname=="Rayl" )  ProcessBit = 3;     
      if( processname=="phot" )  ProcessBit = 4;     
      if( processname=="conv" )  ProcessBit = 5;     
      if(  volumename=="pFid" )  isFid = 1;     
      nt[iev]->Fill(igamma-1,istep,isFid,ProcessBit,xpre,ypre,zpre,xpos,ypos,zpos,edep,ekin);

      istep++;
    }

    if(IsFid) cFid++;

    //-----
    hNPhoton->Fill(nPhoton);
    hNElectron->Fill(nElectron);
    
  }


  fwrite->Write();
  cout<<"nFid = "<<cFid<<endl;

}

void read2(){

  int NO = 4;

  TFile *fread=new TFile("./out_30000.root");

  TGraph *GammaXY=new TGraph();;
  TGraph *GammaXZ=new TGraph();;
  TGraph *GammaYZ=new TGraph();
  TGraph *GammaXYedep=new TGraph();;
  TGraph *GammaXZedep=new TGraph();;
  TGraph *GammaYZedep=new TGraph();

  TPaveText *pave1[30];
  TPaveText *pave2[30];
  TPaveText *pave3[30];
  TPaveText *pave=new TPaveText(0.1,0.1,0.9,0.9);;
  int k=0;    
  int l=0;
  char ProcessName[6][32]={"Other","Transport","Compton","Rayleigh","PhotEleEff","Conversion"};

  int nev=10000;

  for(int i=0;i<nev;i++){

    TNtuple *nt=(TNtuple*)gROOT->FindObject( Form("nt_%d",i) );
    float itrack,istep,isFid,ProcessBit;
    float xpre,ypre,zpre,xx,yy,zz,edep,ekin;
    float ekinpre=0;
    float ekinsum=0;
    int trackpre=0;
    float alledep=0;    

    nt->SetBranchAddress("itrack",&itrack);
    nt->SetBranchAddress("istep",&istep);
    nt->SetBranchAddress("isFid",&isFid);
    nt->SetBranchAddress("ProcessBit",&ProcessBit);
    nt->SetBranchAddress("xpre",&xpre);
    nt->SetBranchAddress("ypre",&ypre);
    nt->SetBranchAddress("zpre",&zpre);
    nt->SetBranchAddress("xx",&xx);
    nt->SetBranchAddress("yy",&yy);
    nt->SetBranchAddress("zz",&zz);
    nt->SetBranchAddress("edep",&edep);
    nt->SetBranchAddress("ekin",&ekin);
    
    int nstep=nt->GetEntries();
    //cout<<nstep<<endl;
    for(int j=0;j<nstep;j++){
      
      nt->GetEntry( j );
      
      if(i==NO){
	  GammaXY->SetPoint(j,xx,yy);
	  GammaXZ->SetPoint(j,xx,zz);
	  GammaYZ->SetPoint(j,yy,zz);
	  if( (int)ProcessBit!=1 ){
	    GammaXYedep->SetPoint(k,xx,yy);
	    GammaXZedep->SetPoint(k,xx,zz);
	    GammaYZedep->SetPoint(k,yy,zz);
	    k++;
	    //	    cout<<ProcessName[(int)ProcessBit]<<" "<<itrack<<endl;
	    
	    if( (int)isFid==1 ){
	      pave1[l]=new TPaveText(xx-20,yy+10,xx+20,yy+30,"user");
	      pave1[l]->AddText( Form("%d",j) );
	      pave2[l]=new TPaveText(xx-20,zz+10,xx+20,zz+30,"user");
	      pave2[l]->AddText( Form("%d",j) );
	      pave3[l]=new TPaveText(yy-20,zz+10,yy+20,zz+30,"user");
	      pave3[l]->AddText( Form("%d",j) );
	      char tmp[64];
	      sprintf(tmp,"No:%d  Process : %s",j,ProcessName[(int)ProcessBit]);
	      pave->AddText( tmp );
	      pave->AddText( Form("X=%.1f  Y=%.1f  Z=%.1f (mm)  Ekine=%.3f (keV)",xx,yy,zz,ekin*1000) );
	      pave->AddText("");
	      //	    cout<< Form("Xpre=%.1f X=%.1f  Ypre=%.1f Y=%.1f  Zpre=%.1f Z=%.1f (mm)  Edep=%.3f (keV)",xpre,xx,ypre,yy,zpre,zz,ekin*1000)<<endl;;
	      l++;
	      
	    }
	  }
	  
	  if(itrack!=trackpre){
	    ekinsum=0;
	    ekinpre=ekin;
	  }
	  else	if(j!=0)  ekinsum += (ekinpre-ekin);
	  //	  cout<<"trackID="<<itrack<<"  "<<ProcessName[(int)ProcessBit]<<"  Ekine="<<ekin<<"  dEkine="<<ekinpre-ekin<<"  dESum="<<ekinsum<<endl;

      }

      ekinpre = ekin;
      trackpre = itrack;
    }
   

    nt->Delete(); 
  }
  

  TCanvas *c1=new TCanvas("c1","c1",1000,800);
  c1->SetGrid();
  c1->Divide(2,2);

  TLegend *leg1=new TLegend(0.6,0.2,0.85,0.35);
  leg1->SetFillStyle(0);

  c1->cd(1);
  TH1 *frame1=gPad->DrawFrame(-180,-180,180,180);
  frame1->SetTitle("XY-Display");
  frame1->SetXTitle("X (mm)");
  frame1->SetYTitle("Y (mm)");
  frame1->Draw();
  c1->cd(2);
  TH1 *frame2=gPad->DrawFrame(-180,-650,180,-300);
  frame2->SetTitle("XZ-Display");
  frame2->SetXTitle("X (mm)");
  frame2->SetYTitle("Z (mm)");
  frame2->Draw();
  c1->cd(3);
  TH1 *frame3=gPad->DrawFrame(-180,-650,180,-300);
  frame3->SetTitle("YZ-Display");
  frame3->SetXTitle("Y (mm)");
  frame3->SetYTitle("Z (mm)");
  frame3->Draw();

  c1->cd(1);
  GammaXY->SetLineColor(4);
  GammaXY->SetMarkerColor(4);
  GammaXY->SetMarkerStyle(21);
  GammaXY->SetMarkerSize(1);
  GammaXY->Draw("pl");
  if(k>0){
    GammaXYedep->SetLineColor(95);
    GammaXYedep->SetMarkerColor(95);
    GammaXYedep->SetMarkerStyle(21);
    GammaXYedep->SetMarkerSize(1);
    GammaXYedep->Draw("p");
  }
  for(int ii=0;ii<l;ii++){
    pave1[ii]->SetFillStyle(0);
    pave1[ii]->SetBorderSize(0);
    pave1[ii]->Draw();
  }

  c1->cd(2);
  GammaXZ->SetLineColor(4);
  GammaXZ->SetMarkerColor(4);
  GammaXZ->SetMarkerStyle(21);
  GammaXZ->SetMarkerSize(1);
  GammaXZ->Draw("pl");
  if(k>0){
    GammaXZedep->SetLineColor(95);
    GammaXZedep->SetMarkerColor(95);
    GammaXZedep->SetMarkerStyle(21);
    GammaXZedep->SetMarkerSize(1);
    GammaXZedep->Draw("p");
  }
  for(int ii=0;ii<l;ii++){
    pave2[ii]->SetFillStyle(0);
    pave2[ii]->SetBorderSize(0);
    pave2[ii]->Draw();
  }

  c1->cd(3);
  GammaYZ->SetLineColor(4);
  GammaYZ->SetMarkerColor(4);
  GammaYZ->SetMarkerStyle(21);
  GammaYZ->SetMarkerSize(1);
  GammaYZ->Draw("pl");
  if(k>0){
    GammaYZedep->SetLineColor(95);
    GammaYZedep->SetMarkerColor(95);
    GammaYZedep->SetMarkerStyle(21);
    GammaYZedep->SetMarkerSize(1);
    GammaYZedep->Draw("p");
  }
  for(int ii=0;ii<l;ii++){
    pave3[ii]->SetFillStyle(0);
    pave3[ii]->SetBorderSize(0);
    pave3[ii]->Draw();
  }

  c1->cd(4);
  pave->SetFillStyle(0);
  pave->SetBorderSize(0);
  pave->SetTextAlign(11);
  pave->Draw();

  TArc *aInner=new TArc(0.,0.,133.,45.,315.);
  TLine *lInner=new TLine(133./sqrt(2),-133./sqrt(2),133./sqrt(2),133./sqrt(2));
  aInner->SetFillStyle(0);
  aInner->SetNoEdges(1);
  aInner->Draw();
  lInner->Draw();
  TLine *lFid1=new TLine(-133,-360,-133,-550);
  TLine *lFid2=new TLine(-133,-550,  94,-550);
  TLine *lFid3=new TLine(  94,-550,  94,-360);
  TLine *lFid4=new TLine(  94,-360,-133,-360);
  TLine *lFid5=new TLine(-133,-550, 133,-550);
  TLine *lFid6=new TLine( 133,-550, 133,-360);
  TLine *lFid7=new TLine( 133,-360,-133,-360);

  c1->cd(1);
  aInner->Draw();
  lInner->Draw();
  c1->cd(2);
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();
  c1->cd(3);
  lFid1->Draw();
  lFid5->Draw();
  lFid6->Draw();
  lFid7->Draw();

  leg1->Draw();



}
